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Abstract 

Guided by molecular dynamics simulations, we generalize the Navier-Stokes-Fourier constitu- 
tive equations and the continuum motion equations to include both transverse and longitudinal 
temperatures. To do so we partition the contributions of the heat transfer, the work done, and 
the heat flux vector between the longitudinal and transverse temperatures. With shockwave 
boundary conditions time-dependent solutions of these equations converge to give stationary 
shockwave profiles. The profiles include anisotropic temperature and can be fitted to molecular 
dynamics results, demonstrating the utility and simplicity of a two-temperature description of 
far-from-equilibrium states. 
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Figure 1: Schematic stationary shockwave. Cold fluid enters at the left cold boundary, with 
speed u s ; hot fluid leaves at the right hot boundary, with speed u s — u p . We choose a coordinate 
frame which moves leftward, at speed u s relative to the laboratory frame. The shockwave 
remains stationary in this coordinate frame. 




-7.0 < x < +7.0 

Figure 2: Stationary shockwave. Snapshot from a 10-row molecular dynamics simulation with a 
periodic height of 10-^/3/4. The simulations analyzed in the text are based on 80-row molecular 
dynamics with a periodic height of 80^3/4. 

I. INTRODUCTION 

A leftmoving piston, impacting a fluid with velocity — u p , generates a leftmoving shock- 
wave with velocity — u s . Throughout this paper we analyze such a shockwave from the 
viewpoint of a coordinate system moving leftward, so as to keep pace with the shock. See 
Figures 1 and 2. In this special uniformly-translating coordinate frame the shockwave is 
stationary, simplifying theoretical analyses. One-dimensional stationary Shockwaves^— 
provide a useful computational laboratory for the study of stationary far-from-equilibrium 
states. In such a shockwave a cold fluid is converted irreversibly to a hot one. As the fluid 
moves from left to right, in the shock-centered coordinate frame of the Figures, at speed 
u(x), the x coordinate increases; typically, the corresponding density, the longitudinal 
component of the pressure tensor, and the energy all increase too, in just such a way that 
the spatial structure of the wave is stationary: 

{ u = x,p,P xx ,e } > , 



(du/dt) x = ; (dp/dt) x = ; (dP xx /dt) x = ; (dP yy /dt) x = ; (de/dt) x = . 
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As the velocity decreases from its leftmost entrance value, u(x — > — oo) = u s , to its 
rightmost exit value, u(x — > +00) = u s — u p , the stationary nature of the wave requires 
that the fluxes of mass, momentum, and energy remain constant throughout: 

(pu) x = (pn) co id = (pu)hot ; 

(Pxx + PU 2 )x = (P + pU 2 )cold = (P + PU 2 )hot ; 

(pu)[{e + (P xx /p) + {u 2 /2)] x + Q X = 

{pu)[e + (P/p) + (M 2 /2)] cold = {pu)[e + (P/p) + (u 2 /2)] hot . 

The notation here is conventional, with the pressure tensor P and heat flux vector Q 
assumed to be calculable from the density p, velocity u, energy e, and their gradients. 

Temperatureii^ 1 ^— is our special interest in this work. Temperature is most simply 
and usefully defined as a velocity fluctuation, the "kinetic temperature" : 

kT xx = m((x - (A)) 2 ) ; kT yy = m((y - (y)) 2 ) . 

The angular brackets imply a local average. The velocities here are individual particle 
velocities, whose local average would be the hydrodynamic flow velocity u. Temperature 
is just the fluctuation about this average. It is evident that T xx and T yy can differ. In 
dilute-gas kinetic theory, the difference corresponds to a shear stress: 

pk(T xx - T yy )/ (2m) = (P xx - P yy )/2 [Dilute Gas] , 

where k is Boltzmann's constant and m is the particle mass, which we choose equal to 
unity in what follows. In dense fluids there is no simple relationship between the two 
tensors so that special evolution equations for T xx and T yy need to be developed, as we 
do in Section III. 

The cold fluid, initially moving to the right at the entrance velocity, or "shock velocity" 
u s , is slowed by its encounter with the wave until it reaches its exit velocity u s — u p , where 
Up is the "piston velocity" or "particle velocity". In this irreversible deceleration the 
kinetic energy lost by the decelerating fluid is converted into additional hot fluid enthalpy 
(H = E + PV h = e + Pv): 

Kot - frcoid = [e + (P/p)] hot - [e + (P/p)]coid = [pco\dU 2 /2] - [phot(« s - u p ) 2 /2] . 
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-3.5 < x < +3.5 -3.5 < x < +3.5 

Figure 3: A snapshot spatial profile of a nominally steady one-dimensional shockwave from 

molecular dynamics, using a short-ranged repulsive potential. Spatial one-dimensional averages 
of the temperatures and heat flux (left) and the pressures, density, and energy (right) have 
been computed with Lucy's weight function using a range h = 3. The cold zero-pressure, zero- 
temperature triangular lattice is compressed to twice the initial density (y/4/3 — > 1\j 4/3) by 
the shockwave, just as in Figure 2. 

The cold and hot boundary conditions enclosing the shock are equilibrium ones imposed 
far from the shockfront so that the small-system surface effects complicating the number- 
dependence of nonequilibrium systems are minimized. In implementing these ideas no 
arbitrary or artificial assumptions have to be made. All the observed phenomena follow 
from the assumed form for the interparticle forces. Figures 3, 4, and 5 show typical results 
from molecular dynamics, as is described in more detail in Section II. Notice that the rise 
in longitudinal temperature T xx can be much larger and can occur somewhat earlier— 
than that of the transverse temperature T yy . 

In Section III we discuss the continuum mechanics of the same shockwave problem. 
Evidently any continuum formulation must first of all include the continuum conservation 
laws for mass, momentum, and energy: 

p = -pV ■ u ; 

pu = -V ■ P ; 

pe = -Vtt : P - V • Q . 

Here the pressure tensor P and heat flux vector Q measure the momentum and energy 
fluxes in the local "comoving" (or "Lagrangian" ) coordinate frame moving with the mean 
velocity u(x). Now the superior dot notation is used to indicate the time derivatives of 
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0.4 < (V/N) < 0.9 0.4 < (V/N) < 0.9 

Figure 4: Volume dependence of the temperature tensor (left) and the pressure tensor (right) in 

the stationary shockwave of Figure 3, as calculated with molecular dynamics. Spatial averages 
have been computed with Lucy's weight function using a range h = 3, as is discussed in Section 
II. 
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-3 < time < +3 -3 < time < +3 

Figure 5: Stationary temporal profile for the one-dimensional shockwave of Figure 3, using a 

short-ranged repulsive potential. Spatial averages of the temperatures (left) and the pressures, 
density, and energy (right) have been computed with Lucy's weight function using a range 
h = 3. The initial stress- free cold triangular lattice is compressed to twice the initial density 
by the shockwave, as in Figure 2. The time origin has been chosen, arbitrarily, close to the 
shockfront. 

p, u, and e following the motion at velocity u. In the continuum description these field 
variables are continuous differentiable functions of space and time so that the spatial 
averaging (necessary to an analysis of molecular dynamics data) is unnecessary. 

The steady nature of the shock process makes it possible to use either space or time 
as an independent variable. On the average, the progress of a particle traveling through 
the shockwave follows from the integral of the flow velocity. To illustrate, consider again 
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Figure 6: Stationary spatial profile for a one-dimensional shockwave according to the usual 
Navier-Stokes-Fourier equations for the model fluid: P cq = pe ; e = (p/2) + kT with unit shear 
viscosity, zero bulk viscosity, and unit Fourier heat conductivity. Here the temperature T (left) 
is a scalar, as in conventional continuum mechanics. 

the molecular dynamics profiles shown in Figures 3, with space as the abscissa. Exactly 
the same profiles can alternatively be expressed with time as the abscissa, as in Figure 5. 
To change from space-based to time-based profiles requires use of the ratio (dx/dt) = u: 



where u(x) is the hydrodynamic flow velocity. Thus all the spatial snapshots or equivalent 
temporal wave profiles catalog the sequence of time-ordered states through which the 
particles in a typical volume (initially at x ) pass as they transit the shockwave. 

Because the conventional Navier-Stokes-Fourier approach, illustrated in Figure 6, as- 
sumes a scalar temperature, T = T xx = T yy , several modifications of the continuum de- 
scription need to be made to model the two-temperature results of Figures 3-5 found with 
molecular dynamics, with T xx ^ T yy . In Section III we describe simple modifications of 
the Navier-Stokes-Fourier constitutive and flow equations, along with a numerical method 
which converges nicely to give stationary shockwave profiles in the two-temperature case. 

Section IV is reserved for a summary and our concluding remarks, including sugges- 
tions for adapting our ideas to detailed two- and three-dimensional descriptions of the 
fluctuations in nonequilibrium systems. 
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II. RESULTS FROM MOLECULAR DYNAMICS 



The molecular dynamics simulations leading to our current results are all based on a 
very simple model two-dimensional system of unit-mass unit-radius particles interacting 
pairwise with a short-ranged normalized repulsive potentiat^ 1 ^: 



The length and energy scales set by this potential correspond to the range and strength of 
the interparticle pair forces. The equilibrium properties for this potential can be approx- 
imated very roughly by a theoretical model (based on a random distribution of particles 
in space) resembling van der Waals' mean-field idea, 



P, p, e, and T are the pressure, density, energy, and temperature. Though the models 
and language here all refer to systems in two space dimensions the same ideas can be 
applied equally well to three-dimensional systems. 

We expect that the nonequilibrium properties for this model will likewise provide a 
simple interpretation. We are particularly interested here in generalizing the notion of 
temperature to the tensor case, T xx ^ T yy . The need for this generalization stems from 
the molecular dynamics shockwave simulations summarized in Figures 3, 4, and 5. 

Stationary Shockwaves were obtained from molecular dynamics by matching the mass 
flux of a cold stress-free lattice (p = a/4/3 and speed 1.930) to the mass flux of the hot 
fluid exiting at the righthand boundary (with p = 2^/4/3 and speed 0.965): 



With this choice for the shockwave speed u s = 1.93 and particle (or piston) speed u v = 
u s /2 the shockwave is stationary and corresponds to twofold compression, a "strong" 
shockwave^. The Mach number M = u/c s is not a useful description here as the sound 
speed c s vanishes in the cold state. The momentum and energy fluxes throughout the 
wave are equal to those of the initial cold lattice: 




P = pe; e = (p/2) + kT . 



pu 



Pcold^cold 



PhotWhot = 1-93 x a/473 = 2.229 . 



P xx + pu 2 = v / 473(l-93) 2 = 4.301 ; 



pu[e + (P xx /p) + (u 2 /2)} + Q X = v / 473(1.93) 3 /2 = 4.151 . 
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Spatial averages within the shockwave were calculated here using Lucy's weight 
function^ 3 -^^, 

w Lucy (\x\ <h) = (5/4/i)[l - 6r 2 + 8r 3 - 3r 4 ] ; r = \x\/h < 1 , 

with a range equal to three times the range of the potential, h = 3. The internal energy 
at a gridpoint coordinate x, for example, is computed as a ratio of sums: 

_ J2jW(x - Xj)ej 
e{x) — , 
2^w(x - Xi) 

where the energy of Particle % is the sum of its kinetic energy relative to the local flow 
velocity u(x) plus half its summed-up interaction energy with other nearby Particles {j}. 

Consider now the results shown in Figures 3 and 4. The density, energy, and pressure 
agree roughly with the hyperbolic-tangent profiles derived by Landau and Lifshitz for 
a weak shockwave with constant transport coefficients 3 . Figure 4 shows the pressure- 
temperature-volume states through which the moving fluid travels. The Rayleigh Line, a 
straightline relation linking P xx and the volume, is necessarily satisfied and corresponds to 
the conservation of momentum. In marked contrast, the molecular dynamics temperature 
shows a strong maximum (as might be expected from the mixing of cold and hot Gaussian 
distributions suggested by Mott-Smith 1 ) at the shockfront. Because the work done in 
compressing the fluid appears first in the longitudinal direction we expect that the rise 
in T xx precedes that of T yy , as is confirmed in Figure 3. This thermal anisotropicity 
differs from the conventional textbook result and is the main motivation for our work on 
a two-temperature continuum description, detailed in the following Section. 

III. RESULTS FROM CONTINUUM MECHANICS 
A. General Considerations 

Continuum models combine the universal conservation laws (mass, momentum, and 
energy) and the corresponding evolution equations (continuity, motion, and energy) with 
specific constitutive models. The constitutive models describe the pressure tensor and 
the heat flux vector for nonequilibrium systems. The usual Navier-Stokes assumptions, 
which we follow here for a two-dimensional fluid, are that the pressure tensor and heat 
flux vector respond linearly to velocity and temperature gradients: 

p = p«i - A[V ■ u]I - rj[Vu + W] ; A = 7] V - V ■ 
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Q = -kVT . 

It needs to be emphasized that the choice of particular expansion variables, here V« 
and VT, affects the solutions of nonlinear problems like shockwave structure. Garcfa- 
Colfn and Green emphasized that the description of nonequilibrium continuum mechanics 
is ambiguous whenever the choice of "equilibrium" variables - energy or longitudinal 
temperature or transverse temperature in this case - is ambiguous 17 . The numerical 
value of a Taylor's series in the deviations from equilibrium, truncated after the first 
nonlinear term, is clearly sensitive to the choice of independent variable. 

In the nonequilibrium pressure tensor the superscript * indicates the transposed tensor 
and / is the unit tensor: 

-^11 — I22 — 1 j I \2 — I21 = ) 

rj is the shear viscosity, and A = rj v — rj is defined by the bulk viscosity rj v . In the shockwave 
problem the pressure-tensor definitions give 

Pxx = P eq ~ (Vv + v)du/dx ; P yy = P cq - (r] v - r])du/dx . 

For a two-temperature continuum model it is necessary to formulate the "equilibrium 
pressure" P eq as a function of the (nonequilibrium) energy, density, and the two temper- 
atures. The viscosities and conductivity could likewise depend upon these state variables 
and k can be a tensor, as we show later, with an example. 

When we define T xx and T yy as continuum state variables it becomes necessary for us 
to formulate constitutive relations for their evolution. The simplest such models begin by 
separating the energy into two parts: a density-dependent "cold curve" e cold (p) and an 
additional kinetic or "thermal" part, proportional to temperature: 

e = e cold (p) + e*"^, T yy ) = e cold + (ck)(T xx + T yy ) , 

where ck is a scalar heat capacity. The functional form of the cold curve produces a 
corresponding contribution to the pressure: 

P cold = -de cold /d{V/N) = p 2 de cold /dp . 

Griineisen's 7 defines a corresponding thermal pressure: 

P thermal thermal 
1 
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The viscous part of the pressure tensor is Newtonian: 

pviscous = _ AV . uI _ r]{ y u + Vu t^ _ 

The thermal and viscous parts of the First-Law energy change are then apportioned 
between the x and y directions so as to be consistent with overall energy conservation: 

e = e e (p) = ckT xx + ckT yy ; 

pckf xx = -aWu : (P - JP cold ) -0V.Q + pck(T yy - T xx )/r ; 

pcA;T ra = (a - 1)V« : (P - /P cold ) + (/3 - 1)V • Q + - T yy )/r . 

The thermal relaxation time r has been introduced in the evolution equations to guarantee 
thermal equilibrium far from the shockwave: 

K x ^ ^ ^xx ^yy ^ 

In what follows we consider two models for the cold curve and the heat capacity. First, 
a weak repulsive pair force suggests implementing a "van der Waals model" : 

e cold = (p/2) ; e thcrmal = k(T xx + T yy )/2 ; P c * = pe 

Second, a triangular-lattice-based model, based on Griineisen's ideas, uses the nearest- 
neighbor static lattice energy and pressure corresponding to the pair potential evaluated 
at the nearest-neighbor lattice spacing r, = (10/7r)(l — r) 3 : 

e cold = (30/tt)(1 - r) 3 ; p cold (V/N) = (45/7r)r(l - rf ; 



The corresponding equilibrium equation of state separates the energy and pressure into 
"cold" and "thermal" parts: 

^eq g co '^ -)- gthcrmal . peq pcold _|_ ^gthcrmal 

with 7 chosen so as to roughly reproduce equation of state data from molecular dynamics. 
Let us next apply these two simple cold-curve models to the shockwave problem. 
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B. Potential plus Kinetic van der Waals Models 



First consider an arbitrary, but simple and natural, choice: 

p«i = pe ; e eq = e cold + e thcrmal = (p + kT xx + kT yy )/2 . 
pcold = pe coid = p 2j 2 ^ 

with an initial density of unity and an initial temperature of zero. Twofold compression 
of the cold van der Waals fluid gives the following solution relating the initial and final 
equilibrium states: 

p : 1 ->■ 2 ; u : 2 — y 1 ; T : ->■ 1/4 ; e : 1/2 ->■ 5/4 ; P : 1/2 ->■ 5/2 . 

The mass, momentum, and energy fluxes connecting these states must be constant 
throughout the profile: 

pu = 2 ; P xx + pu 2 = 9/2 ; pu[e + (P xx /p) + (u 2 /2)} + Q x = 6 . 

Consider the most extreme anisotropic situation consistent with energy conservation, 
in which all the work done and heat transfered are associated with thermal change in the 
x direction. The thermal relaxation time r, here chosen equal to unity, guarantees that 
the x and y temperatures equilibrate in a time of order r: 

•thermal = • _ . C old (p) = {k/2 )(f xx + f yy ) ; 

p(k/2)f xx = -Vu : (P - JP cold ) - V • Q + p{k/2){T yy - T xx )/t ; 

p(A;/2)T OT = p(k/2)(T xx - T ro )/r ; r = 1 . 

Two solutions of these equations appear in Figures 7 and 8. For both of them we chose 
a shear viscosity of unity and a vanishing bulk viscosity: 

Pxx = P cq - du/dx ; P yy = P cq + du/dx . 

The heat flux vector requires that an additional choice be made for its response to 
the gradients of T xx and T yy . We compare two choices in Figures 7 and 8. For both of 
them the overall conductivity is unity, but the heat flux responds differently to the two 
components of VT: 

Q x = —nVTyy = —VT yy [Choice 1] . 
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-5.0 < x < +5.0 -5.0 < x < +5.0 

Figure 7: Typical solution of the generalized Navier-Stokes-Fourier equations for the van der 

Waals model with both heat and work contributing to T xx and with the heat flux responding 
only to the gradient of T yy . The shear viscosity, heat conductivity, heat capacity, and thermal 
relaxation times are all taken equal to unity. 

Q x = ~k(VT xx + VT yy )/2 = -(VT XX + VT yy )/2 [Choice 2] ; 

It is good fortune that the Shockwave equations we summarize here are relatively easy 
to solve numerically. The usual numerical method is the "backward Euler" scheme^. One 
starts near the "hot" boundary and integrates backward, using a first-order difference 
scheme. That approach fails here, due to the temperature relaxation terms, which are ex- 
ponentially unstable in the time-reversed case. An integration forward in time is required 
in the presence of relaxation. A successful "staggered-grid" (two separate spatial grids) 
algorithm results if the density p c is defined at cell centers and energy, temperature, and 
pressure are defined at the nodes which bound the cells^ 1 ^. This algorithm follows the 
dynamics correctly and converges nicely to the stationary profiles shown in Figures 7 and 
8. A computational mesh spacing of dx — 0.1 is sufficient, using the second-order spatial 
differencing scheme outlined in References 18 and 19 with fourth-order Runge-Kutta time 
integration. 

In the early days of shockwave modeling this computational simplicity was by no means 
apparent, so that there is an abundant literature on the stability of numerical methods 
for the shockwave problem^. Now, in the early days of tensor-temperature models, the 
main challenge is to develop well-posed constitutive equations consistent with both the 
conservation laws and the empirical results from molecular dynamics. 

Interesting aspects of both solutions are (i) the minimum in P yy (x), which suggests the 
need for bulk viscosity in modeling molecular dynamics results, and (ii) the pronounced 
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-5.0 < x < +5.0 -5.0 < x < +5.0 

Figure 8: Typical solution of the generalized Navier-Stokes-Fourier equations for the van der 

Waals model with both heat and work contributing to T xx and with the heat flux responding 
equally to the gradients of both T xx and T yy . The shear viscosity, heat conductivity, heat 
capacity, and thermal relaxation times are all taken equal to unity. 

maximum in T xx (x), leading the response of T yy and roughly equal in magnitude to that 
found in the dynamical results of Section II. 

The physical ideas incorporated in this simplest approach are four: (i) the pressure 
and the work done can usefully be separated into a "cold" part and a "thermal" part; 
(ii) the heat flux Q responds to a linear combination of the temperature gradients VT XX 
and VTyy in the usual way, supplemented by (iii) the thermal relaxation of the thermal 
anisotropicity, and (iv) separate linear combinations of the work done and heat absorbed 
contribute to T xx and T yy throughout the shock compression process. 

Here the total pressure, P = P $ + P , contains potential and kinetic components, 
measurable separately with molecular dynamics. These extensions of the Navier-Stokes 
approach closely parallel the relaxation-time treatments of strong ideal-gas Shockwaves 
carried out by Xu, Josyula, Holian, and Mareschal^^. Our more general approach 
necessarily differs from theirs by allowing for contributions from the potential energy to 
temperature changes and the transfer of heat. The pressure profiles shown in Figures 7 
and 8 also indicate the need for bulk viscosity, in that the molecular dynamics results show 
a monotone- increasing P yy , in contrast to the distinct minimum found here in the absence 
of bulk viscosity. We turn next to a slightly more sophisticated model, an extension of 
Griineisen's equilibrium equation of state. 
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C. Cold plus Thermal Griineisen Models 



For gases, where the pressure and temperature tensors are proportional to one another, 
a systematic expansion of the Boltzmann equation can be, and has been, triecU&ii^ 1 ^. 
Xu and Josyulaii as well as Holian and Mareschal^ developed solutions of generalized 
relaxation-time Boltzmann equations for the shockwave problem. For dense fluids only 
Enskog's hard-sphere-based theory is available. More flexible empirical models need to be 
developed for dense-fluid Shockwaves. A trial set of two-temperature evolution equations, 
the simplest plausible set generalizing the van der Waals model above, makes use of 
Griineisen's "cold curve" representation of the energy and pressure to define "thermal" 
contributions. These thermal parts include both the effects of thermal agitation (heat 
and temperature) and of mechanical distortion (work, through compression with viscous 
deformation) : 

Z7> /f-cold i rpthermal tj 7->cold , othermal , rjviscous 

£j = <P + -C/ , f[x X an( j yy ) = f + f + f 

For the molecular dynamics simulations discussed in Section II the cold parts of the 
pressure and energy, as well as their time dependence, are naturally defined by imagining 
a perfect static triangular lattice of particles: 

E cold /N = e cold = (30/tt)(1 - r) 3 ; P cold V/N = -(dE cold /dV) = (45/vr)r(l - r) 2 . 

p gCold = . pcold _ 

Here r is the separation of the six nearest neighbors in a cold triangular lattice, so that 
p = 

Just as in the equilibrium Griineisen model the thermal energy and the nonviscous 
part of the thermal pressure are taken to be proportional to temperature: 

e thermal = + Kj/flf ; pthermal = ^thermal ? 

where 7 is Griineisen's constant and ck is a heat capacity. 

The Krook-Boltzmann relaxation terms, with relaxation time r, are the simplest means 
for guaranteeing thermal equilibrium, with the two temperatures approaching one another 
far from the shockfront. 

Because molecular dynamics simulations indicate that temperature becomes a tensor in 
strong shockwaves, a tentative two-temperature formulation can be based on separating 
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the internal energy and the pressure into the three components suggested by classical 
statistical mechanics, including Newtonian shear and bulk viscosities: 

E = Ne = $ cold + $ thcrmal + K x + K y ; 

Pxx = P eq - (v + r] V )du/dx ; P yy = P eq + (rj - r) v )du/dx ; 
P eq = p[0 cold + 1C kT xx ] or p[0 cold + 1C kT yy } or p[0 cold + 1C k(T xx + T yy )/2\ ; 

formal = ^thermal + ( k/2 )(T xx + T yy ) = ck(T xx + T yy ) . 

The sum of the three energy evolution equations just given is designed to reproduce 
the usual First Law energy equation, 

E = Eq — E w , 

where Eq and E w are the comoving rates at which heat enters the fluid and at which the 
fluid performs work on its surroundings. The constitutive relations for P and Q must also 
be given. For a two-dimensional Newtonian fluid with shear viscosity rj and bulk viscosity 
i] v we have 

Pxx = P cq - (v + Vv)du/dx ; P yy = P cq + (rj - r) v )du/dx . 

The heat flux is given by a generalization of Fourier's law, with independent contributions 
from VT XX and VT yy . 

Additional generalizations of this approach can be developed as needed to describe re- 
sults from simulations. It is only required that any such model satisfy energy conservation 
and reduce to the Navier-Stokes- Fourier model in the weak-shock limit. To illustrate the 
possibilities, compare the molecular dynamics results of Figure 3 to the model calculations 
of Figure 9. In Figure 9 the relaxation time has been increased to 3, the heat capacity 
doubled, to ck = 2k, and the heat conductivity set equal to 6 so as to better match the 
empirical results of molecular dynamics. The value of Gruneisen's 7 is 0.3, and the bulk 
and shear viscosities are both equal to unity. The results from these choices (which are by 
no means optimized) resemble the shockwave profiles obtained with molecular dynamics. 

IV. CONCLUSIONS AND PROBLEMS FOR THE FUTURE 



We have shown here that it is relatively easy to model the thermal anisotropicity found 
in atomistic simulations of strong Shockwaves. Thermal relaxation, bulk viscosity, and 
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Figure 9: Solution of the generalized Navier-Stokes-Fourier equations with both heat and work 

contributing solely to T xx and with the heat flux Q = —K(VT xx -\-TVTyy)/8. The shear viscosity, 
bulk viscosity, heat conductivity, and thermal relaxation times are respectively 1, 1, 6, and 3. 
Griineisen's 7 is 0.3 and ck = 2k. 

Griineisen equations of state are useful components of a kinetic shockwave model. By 
apportioning the longitudinal and transverse thermal portions of the work, heat, and 
heat flux vector a variety of useful models can be developed and used to reproduce results 
from simulations. A forward-in-time fourth-order Runge-Kutta (as opposed to backward 
Euler) integration of the cell and nodal motion equations results in accurate and stable 
continuum dynamics. 

One of the recent observations from molecular dynamics is that the stress and heat 
flux lag somewhat behind the strainrate and the temperature gradient^ 2 . It is desirable 
that models be generalized to reflect these lags. Some study of time-delayed differential 
equations is necessary to model this phenomenon. 

A significant goal is the extension of these same ideas to the fluctuating stress and heat 
flows of two and three dimensional fluids. A comparison of results from molecular dynam- 
ics with those from two and three-dimensional two-temperature continuum simulations 
should provide useful tools for describing fluctuations within the overall one-dimensional 
flows. 

These results show that even far-from-equilibrium shocks can be treated in a semi- 
quantitative way by relating the tensor parts of the energy flows to one another in a 
relatively simple way. An intriguing result of some model calculations is the stable rever- 
sal of the direction of the heat flux vector. Though this reversal seems unphysical, there 
is no difficulty in obtaining stable numerical profiles which include flux reversal. 
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Figure 10: Two snapshots of the collision of two 1600-particle slabs (periodic in the y direction, 
with height 20 and initial width 20^3/4. The initial velocities, u p = ±0.965 give twofold shock 
compression, followed by a nearly isentropic free expansion at the free surfaces. 

The thermodynamic irreversibility of the shockwave process has an interest indepen- 
dent of the definition of temperature and is worth futher study. The shock process it- 
self obeys purely Hamiltonian mechanics, and Liouville's Theorem^. Even so, by using 
Levesque and Verlet's integer version of the leapfrog algorithm^ the entire shockwave dy- 
namics can be precisely reversed, to the very last bit. The apparent paradox, a perfectly 
time-reversible but thermodynamically irreversible process, can most clearly be illustrated 
by simulating the (inelastic) collision of two zero-pressure blocks of fluid. The collision 
of the blocks, with velocities ±u p generates two Shockwaves, with velocities ±(u s — u p ). 
Two snapshots from such a simulation are shown in Figure 10. 
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